SUBROUTINE ZMUSIC_SPECTRUM(SAMPLES, LSAMPLES, N, M, PSPECTRUM, LPSPECTRUM)
  IMPLICIT NONE
  INTEGER :: LSAMPLES, N, M, I, K, LPSPECTRUM, INFO
  COMPLEX*16 :: SAMPLES(LSAMPLES), R(M,M), G(M, M-N), A(M), RES(M)
  DOUBLE PRECISION :: PSPECTRUM(LPSPECTRUM)
  DOUBLE PRECISION :: PI, OMEGA, DELTA, NORM, DZNRM2

  CALL ZCORREST(SAMPLES, LSAMPLES, M, R)

  CALL ZMUSIC_COMMON(R, N, M, G)

  PI = 3.1415926535897932384626433832795028841971d0
  OMEGA = -PI
  DELTA = 2*PI/(LPSPECTRUM-1)
  DO I = 1, LPSPECTRUM
    DO K = 0, M-1
      A(K+1) = EXP((0, -1) * OMEGA  * K)
    END DO
    CALL ZGEMV('C', M, M-N, 1.0d0, G, M , A, 1, 0, RES, 1)
    NORM = DZNRM2(M, RES, 1)
    PSPECTRUM(I) = 1.0d0 / NORM
    OMEGA = OMEGA + DELTA
  END DO
END SUBROUTINE ZMUSIC_SPECTRUM
